Bulges or Bars from Secular Evolution? 
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ABSTRACT 

We use high resolution colhsionless A^-body simulations to study the secu- 
lar evolution of disk galaxies and in particular the final properties of disks that 
suffer a bar and perhaps a bar-buckling instability. Although we find that bars 
are not destroyed by the buckling instability, when we decompose the radial den- 
sity profiles of the secularly-evolved disks into inner Sersic and outer exponential 
components, for favorable viewing angles, the resulting structural parameters, 
scaling relations and global kinematics of the bar components are in good agree- 
ment with those obtained for bulges of late-type galaxies. Round bulges may 
require a different formation channel or dissipational processes. 

Subject headings: galaxies: bulges - galaxies: evolution - galaxies: formation - 
galaxies: kinematics and dynamics - galaxies: photometry - galaxies: spiral 



- 2 - 



1. Introduction 

Evidence has accumulated in the past decade showing that many bulges, especially at 
low-masses, have a disk-like, almost-exponential radial fall-off of the stellar density (An- 
dredakis & Sanders 1994; Courteau et al. 1996; de Jong 1996; CaroUo et al. 1998, 2001; Car- 
oUo 1999; MacArthur et al. 2003), and in some cases disk-like, cold kinematics (Kormendy 
1993; Kormendy et al. 2002). Comparisons of bulge and disk parameters have furthermore 
shown a correlation between the scale-lengths of bulges and disks (de Jong 1996; MacArthur 
et al. 2003) and, on average, similar colors in bulges and inner disks (Terndrup et al. 1994; 
Peletier & Balcells 1996; Courteau et al. 1996). The disk-like properties of bulges and the 
links between bulge and disk properties have been suggested to indicate that bulges may 
form through the evolution of disk dynamical instabilities such as bars, which are present in 
about 70% of nearby disk galaxies (Knapen 1999; Eskridge et al. 2000). 

It has long been known that bars lead to angular momentum redistribution and to an 
associated increase in the central mass density. Already Hohl (1971) found that an initially 
single component disk evolved a double exponential density profile under the infiuence of a 
bar; as a result, the scale- length of the outer disk increases. It has been suggested that bars 
may be efficient at building three-dimensional stellar bulge-like structures via scattering of 
stars at vertical resonances (Combes et al. 1990), or by the coUisionless buckhng instabihty 
(which weakens the bar; Raha et al. 1991), or via bar destruction due to the growth of a 
central mass concentration (Pfenniger & Norman 1990; Norman et al. 1996). As buckling 
occurs relatively easily, it could be a particularly promising avenue for bulge formation. This 
instability leads to the large-scale, coherent bending of the bar perpendicular to the plane of 
the disk (Raha et al. 1991). Buckling, which is a result of vertical anisotropy (Araki 1985; 
Fridman & Polyachenko 1984; Merritt & Hernquist 1991; Merritt & Sellwood 1994), thickens 
the stellar system and weakens the bar. Evidence that buckling occurs in nature has rehed 
on the fact that buckled bars are boxy/peanut-shaped when viewed edge-on and on the 
unmistakable gas-kinematic signature of a bar observed in such bulges (Kuijken & Merrifield 
1995; Merrifield & Kuijken 1999; Bureau & Freeman 1999). The relevance of buckling to 
bulge formation remains however rather anectdotal, as little work has been done to check 
that the structural and kinematic properties of buckled bars are quantitatively consistent 
with those observed in bulges. 

In this Letter, we report on high resolution A^-body simulations of bar-unstable disks, 
some of which buckle, and compare their secularly-evolved structural and kinematic prop- 
erties with those of bulges in local galaxies. We describe our simulations in §2, present our 
results in §3 and discuss their implications in §4. 
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2. Methods 

2.1. Rigid halo simulations 

Most of the A^-body simulations reported in this Letter consist of a live disk inside a 
rigid halo, which allows large numbers of particles and high spatial resolution. The rigid 
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halos were represented by either a logarithmic potential, ^l — ^ l^('^^ + '^h)) ^ Hernquist 
potential $// = — The initially axisymmetric disks were modeled by the generalized 
surface density profile /(r) oc exp[— (r/ro)*^^/")] introduced by Sersic (1968). This gives a 
de Vaucouleurs (1948) profile for n = 4 and an exponential profile for n = 1; we used 
1 < n < 2.5. The disks have scale-length Rd, mass M^, Gaussian thickness and are 
truncated at a radius Rt- Disk kinematics were set up using the epicyclic approximation to 
give constant Toomre-Q, which we varied from 1.2 to 2.0. Vertical equilibrium was obtained 
by integrating the vertical Jeans equation. The disks were represented by 4 — 7.5 x 10^ 
equal-mass particles. We use units where Rd — M — G — 1; thus the unit of time is 
{Rl/GMy/''. 

The simulations were run using a 3-D cyhndrical polar grid code (Sellwood & Valluri 
1997) with Nr X X ^ 60 X 64: X 243. The radial spacing of grid cells increases 
logarithmically from the center, reaching to ~ 2i?t; in all cases, Rt — 5. The vertical 
spacing, 5z, of the grid planes was set such that 0.2 > -2d > 0.025 > iSz. We used Fourier 
terms up to m = 8 in the potential, with a Plummer softening length e = 0.017. We verified 
that our results are not sensitive to resolution by running the most strongly buckled model 
with m = 32 and smaller 6z. Time integration was performed with a leapfrog integrator 
with a fixed time-step, 6t = 0.01 for all runs except for the n > 1 simulations, for which we 
used St — 0.0025. We chose values for the halo parameters such that our rotation curves 
were approximately flat to large radii. For the logarithmic halos, we set Vh — 0.68, while 
the Hernquist halos had Mh = 43.4. We measured the amplitudes of the bar, A^, and of 
buckling, Az, as the normalized amplitudes of the m — 2 tangential and vertical density 
distributions. We will discuss the effects of different initial conditions elsewhere. 



2.2. Modeling and measured parameters 

In order to compare our final systems with observed galaxies, we first obtained, for all 
simulations, radial density profiles of the disk mass density distribution at three inclinations, 
i — (face-on), i — 30° and 60°. For i — 0, we simply computed the azimuthally averaged 
mass profile. For i = 30° and 60°, we considered 3 orientations of the bar with respect to 
the major-axis (0bar — 0,45° and 90°); we then measured the 1-D projected surface density 



-4- 



profiles with the task ellipse in IRAF^. 

We decomposed these mass profiles into a central Sersic component plus an outer expo- 
nential disk; a Sersic law is usually used in modehng bulge fight profiles (Andredakis et al. 
1995; Graham 2001; MacArthur et al. 2003). We stress that by decomposing the final density 
profiles into a "bulge" plus disk component we are not implying that secular evolution has 
produced three-dimensional bulges. Nonetheless, for brevity, we will refer to as "bulges" the 
central Sersic components, and as "disks" the exponential components. 

Our Sersic bulge plus exponential disk decompositions are characterized by five pa- 
rameters: So,(i, So,6, the disk and the bulge central surface densities, respectively; -R^, the 
scale-length of the outer exponential disk; Rb,effj the half-light radius of the bulge; and rif,, 
the index of the Sersic profile. In light of the ill-conditioned fitting described by MacArthur 
et al. (2003), was held fixed with respect to the remaining four parameters when search- 
ing for the best-fitting models; the best-fitting Ub was then identified as the one giving the 
smallest when repeating the fits with 0.1 < < 4 in steps of 0.1. The five parameters de- 
rived from the bulge/disk decompositions allow us to compute three dimensionless structural 
quantities to compare with observations: R^^effl Rd and B/D, the ratio of bulge to disk 
mass (we assumed a constant mass-to-light ratio when comparing with the light-weighted 
measurements available for real galaxies). 

We also measured the bulge eUipticity, 65, the fine-of-sight velocity dispersion, a (by 
averaging the corresponding profiles within some radial range), and 1^, the peak line-of-sight 
velocity within the same radial range on the disk major-axis. These allow us to compute 
Vp/a io investigate the kinematic properties of the resulting bulges in the Vp/a versus eb 
plane. In this plane, normal bulges are thought to follow the locus traced by isotropic oblate 
rotators (Binney & Tremaine 1987), and bulges that result from the secular evolution of 
disks are thought to emerge as systems that are dynamically colder than isotropic oblate 
rotators (Kormendy 1993; Kormendy et al. 2002). 

2.3. Live halo simulation 

Rigid halo simulations are better suited to systems in which the disk is dominant in 
the inner regions (as are the simulations discussed here), because the interaction with the 
halo is then weaker. While massive disks represent a reasonable assumption for high surface 



^IRAF is distributed by NOAO, which is operated by AURA Inc., under contract with the National Science 
Foundation 
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brightness barred galaxies (Debattista & Sellwood 2000), it is important to verify that the 
full interaction with a live dark matter (DM) halo does not lead to a drastically different 
evolution for the structural and kinematic parameters. The ability of a bar to lose angular 
momentum to a DM halo helps make it stronger (Debattista & Sellwood 2000). We therefore 
have run a pair of matched live and rigid halo simulations. The live halo simulations were 
run on PKDGRAV (Stadel (2001); details of the live halo simulations will be presented 
elsewhere) . We checked that these results are not sensitive to resolution by running the hve 
halo model with larger N and smaller e. The two simulations evolved differently: the bar in 
the live halo case was slightly stronger and slows down because of dynamical friction with the 
halo. Nonetheless, the mass density and V^/cr,^ profiles in the rigid- and live-halo simulations 
remain quite similar, as we illustrate in Fig. 1 (we use o"^ in this comparison because it is 
representative of what we might expect in high inclination observations). The larger central 
density in the live halo simulation of Fig. 1 is a result of the fact that its bar can shed 
additional angular momentum. Otherwise, this comparison suggests that the gross features 
seen in the live halo simulation are well-reproduced by the rigid halo simulation, giving 
confidence that the rigid halo simulations survey from which we draw the main conclusions 
of this paper is adequate to describe the evolution of massive disks. Moreover, by using 
rigid halos, we assure that we obtain a minimal level of bar secular evolution, which helps 
to disentangle the effects of different evolutionary processes. 
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Fig. 1. — Comparison of the evolution of the density profiles in matched live (solid lines) 
and rigid halo (dashed lines) simulations after ~ 2.4 Gyr. In the top panel, the initial profile 
is indicated by the dotted hue. The fractional change in the profiles, relative to the initial 
profile, is shown in the middle panel. The bottom panel shows the profile of Vja^. 
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3. Results 

3.1. Time evolution 

The evolution of a representative strongly buckling simulation, shown in Fig. 2, pro- 
duced a bar by i ~ 50, which then buckled strongly at i ~ 100. Despite the strong buckhng, 
the bar is weakened but not destroyed. As is well known, the process of bar formation drives 
a redistribution of angular momentum, leading to an increase in the central density. This 
process is largely complete by the time the bar buckles: the buckling instability does not 
alter significantly the scale length or the mass of the inner Scrsic component. However, it 
is interesting that at buckling, changes from ~ 1 to n?, ~ 1.5: thus buckling may 
contribute to the scatter of around the exponential value observed in the bulges of real 
intermediate-type galaxies. 
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Fig. 2. — The evolution of a representative strongly buckling simulation. Prom bottom to 
top the panels show evolution of (strength of bar; solid hne) and (strength of buckhng; 
dashed line), Rb^eff/Rd, ni, and B/D, all measured a,t i — 0. The bar forms at t ~ 50 and 
buckles at t ~ 100. Before t = 30 the fitting algorithm is unstable and spurious results are 
obtained because the algorithm fits transient spiral structure. 
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Fig. 3. — The structural properties of our simulations compared with observed galaxies 
(shaded area) from MacArthur et al. (2003) and Graham (2003). Squares are filled (empty) 
when ef, < ed (e;, > e^) and are red (green) when buckling is strong (weak). 



3.2. Comparisons with observations 



In Fig. 3 we compare the structural properties of our final systems with similar quantities 
measured for local bulges. Our simulations are represented by filled symbols when is larger 
than the ellipticity of the disk, = 1 — cosi. All i < 30° projections have e;, > e^, as are 
about half of the projections at i = 60°. Bulges in our simulations cover a broad range 
of parameter space, including regions where no real bulges are found. This mismatch is 
diminished when only systems with < are considered, but at the cost that a large fraction 
of the projections, including all at i < 30°, are excluded from the sample. Nevertheless the 
overlap between the simulations and the observations is quite good. 
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Fig. 4. — The V-pjo versus diagnostic plane. Symbols for simulations are as in Fig. 3 while 
the observational data of Kormendy (1993) are indicated by black stars. 
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We also checked whether global kinematics can distinguish between bona-fide bulges 
and bulge-like bars. In Fig. 4 we plot Vp/a versus for real bulges and for our simulations. 
For real bulges, Vp and a were obtained by fitting a de Vaucouleurs profile to the bulge. Our 
simulations arc rather poorly fit by a de Vaucouleurs profile, which introduces a systematic 
difference between our measurements and Kormendy's (1993). To quantify some of this 
uncertainty, we have measured Vp and a within Rb,eff and Rb^efff^ and used half the 
difference as our error estimate. This diagram shows that evolved bars can show, under some 
viewing angles, global kinematic-elongation properties typical not only of the dynamically- 
cold "pseudo-bulges" , but also those that are taken as defining the bona-fide, dynamically- 
hot bulges. These are the expected projection effects for triaxial spheroids seen at high 
inclination (Binney & Tremaine 1987). Thus the population of bulges which arc identified 
by bulge/disk decompositions or by kinematics may have some contribution from bars. 



4. Discussion 

Although Raha et al. (1991) only reported bar-weakening, not destruction, by the buck- 
ling instability, it has become common wisdom that buckling destroys bars. Our simulations 
have failed to turn up a single instance in which the bar was destroyed by buckling, although 
we cannot exclude that it is in some extreme case. The channel of bulge formation by the 
dissipationless destruction of bars during buckling, therefore, is not viable. 

However, the minimal collisionless secular evolution present in our simulations must 
also occur in nature, resulting in systems that exhibit double component mass density pro- 
files. For certain viewing orientations, the spread in structural parameters and kinematic 
properties are indistinguishable from those observed in systems that are classified as bulges. 

Nonetheless, the existence in nature of round bulges inside low inclination galaxies, 
which our simulations cannot reproduce, requires that other processes are also involved. 
Possibly secular evolution including dissipativc gas (Mayer & Wadslcy 2004) will result in 
rounder bulges, but this requires extended central objects with masses of ~ 10 — 20% that 
of the disk (Shen & Sellwood 2004). The higher interaction rate in the early universe may 
have triggered the large gas infiows needed to build such objects. 

Finally, because the halos of our simulations are rigid, the baryonic component cannot 
but conserve its angular momentum. In (semi-) analytic models of disk galaxy formation 
(Fall & Efstathiou 1980; Mo et al. 1998; Lacey & Silk 1991; White & Frenk 1991; Dalcanton 
et al. 1997; van den Bosch 1998), the distribution of disk scale- lengths, Rd, is set by that 
of their angular momenta, de Jong & Lacey (1996) found that the width of the observed 
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Rd distribution at fixed luminosity is smaller than that predicted by analytic theory. Our 
simulations show that, under the influence of a bar, may increase by a factor of 2 or 
more at constant global angular momentum. Since less extended disks arc likely to be more 
bar-unstable, the secular evolution of these disks may be responsible for at least part of this 
discrepancy. 

We would like to thank the anonymous referee for useful suggestions. 
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